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ABSTRACT 

We construct a radiation-hydrodynamics model for the obscuring toroidal 
structure in active galactic nuclei. In this model the obscuration is produced at 
parsec scale by a dense, dusty wind which is supported by infrared radiation pres- 
sure on dust grains. To find the distribution of radiation pressure, we numerically 
solve the 2D radiation transfer problem in a flux limited diffusion approximation. 

We iteratively couple the solution with calculations of stationary ID models 
for the wind, and obtain the z-component of the velocity. 

Our results demonstrate that for AGN luminosities greater than 0.1L e dd 
external illumination can support a geometrically thick obscuration via out- 
flows driven by infrared radiation pressure. The terminal velocity of marginally 
Compton-thin models (0.2 < tt < 0.6). is comparable to or greater than the 
escape velocity. In Compton thick models the maximum value of the vertical 
component of the velocity is lower than the escape velocity, suggesting that a 
significant part of our torus is in the form of failed wind. 

The results demonstrate that obscuration via normal or failed infrared-driven 
winds is a viable option for the AGN torus problem and AGN unification models. 
Such winds can also provide an important channel for AGN feedback. 


1. Introduction 

The active galactic nucleus (AGN) unification scheme envisages the presence of a ge- 
ometrically and optically thick, torus-like structure which wraps and hides a supermassive 
black hole (BH) and active parts of an accretion disk. The paradigm relies on the property of 
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such a structure to obscure the central regions of AGN in type II objects, making the torus 
responsible for the apparent dichotomy of active galaxies (e.g. Antonucci & Miller (1985)). 

Direct evidence for the existence of the toroidal obscuration comes from interferometric 
mid-infrared observations of the nearby Seyfert II galaxies such as NGC 1068, (Jaffe et al. 
2004), and the Cireinus galaxy (Tristram et al. 2007). Studies such as these support the idea 
of a cold (T=100 - 1000 K) torus situated approximately 1 pc away from a supermassive 
BH. These observations also reveal the inner hot (~ 800 K) funnel of the torus and the 
outer, colder (~ 300 K) dusty component (Raban et al. 2009; Bock et al. 2000). Theoretical 
modeling (Krolik & Begelman 1986; Dorodnitsyn et al. 2008) also predicts that the torus 
funnel is significantly hotter than the rest of the torus body due to heating by X-rays 
generated in the inner parts of an accretion disk. 

Indirect evidence for the geometrically thick obscuring structure located at parsec scales 
comes from observations of warm absorber gas. Such observations of nearby Seyfert I galaxies 
by the grating spectrographs on the X-ray telescopes Chandra and XMM-Newton reveal rich 
X-ray line spectra in the 0.1-10 keV range, which contain numerous lines from ions such as 
Fe, Si, S, O, Mg, and Ne, broadened and blue-shifted by 100 - 1000 kms' 1 . These have 
been detected from approximately half of low-redshift AGN (Halpern 1984; Kaspi et al. 2002; 
Steenbrugge et al. 2005; Reynolds 1997; McKernan et al. 2007). Numerical modeling shows 
that if the cold gas of the torus is exposed to extensive X-ray heating then an evaporative 
flow is formed. Simulations suggest that this gas is producing the warm absorber spectrum 
(Dorodnitsyn et al. 2008; Dorodnitsyn & Kallman 2009; Chelouche 2008). 

The wind scatters radiation from the accretion disk and broad-line region toward the 
observer, giving rise to polarized radiation flux observed in the optical and UV (Antonucci 
& Miller 1985). and predicted by theoretical modeling to exist in X-rays in the 0.1 — 10 keV 
range (Dorodnitsyn & Kallman 2010). 

One of the major problems which must be addressed by a theory of AGN obscuration 
is hew the torus resists collapse into a geometrically thin disk. If the torus is supported by 
rotation and gas pressure then the temperature of the gas should be of the order of the virial 
temperature r vir g = 2.6 x 10° Me/r pc K. where M 6 is the BH mass in 10 6 M®, and r pc is the 
distance in parsecs. Clearly, such temperatures cannot be reconciled with the existence of 
dust. 

One solution to this problem is that the pressure of infrared photons on dust prevents 
the vertical collapse of the torus and supports its geometrical thickness. Comparing the 
energy density of the X-ray and UV-photons. 3.44 x 10” 5 M 6 /rp C ergcm~ 3 , (assuming that 
the black hole radiates at half of its Eddington luminosity and a 30% covering fraction of the 
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Compton thick portion of the torus) with the energy density of infrared radiation, we obtain 
that a gas-dust temperature of a fewxlOOK is required if all these X-ray and UV photons 
are converted to the infrared. A more elaborate treatment (see Section 3) shows that if the 
temperature of the torus is a fraction of T V i r , r ~ 527 ( -/IS- ) 1 / 4 _ 938 ( n / 10 ,)t/4 K, where n 
is the number density, the torus thickness will be maintained by radiation pressure. Here 
T v j r , r is another definition of the virial temperature based on the radiation energy density of 
black-body radiation in a radiation-dominated plasma. 

Alternative scenarios assume obscuration either from a warped disk or via a magnetically- 
driven accretion disk wind. The first scenario (Phinney 1989) implies that the transition 
accretion disk is locally geometrically thin but strongly warped (Sanders et al. 1989). In a 
magnetically-driven wind scenario (Konigl & Kartje 1994), the torus is identified with the 
outer regions of a dense hydromagnetic outflow’. 

All these models, including the one w’hich we propose in this paper, describe the torus 
as consisting of tenuous plasma. Regardless of the mechanism for the obscuration, the gas 
of the torus is self-gravitating, susceptible to various instabilities, and so possibly clumped 
or/and in the form of clouds (e.g Elitzur (2008)). 

A global solution requires modeling of the radiatively supported torus via multi-dimensional 
and multi-group radiation hydrodynamics simulations including self-gravitation. To accu- 
rately treat all the macro- and micro-physical processes known to be involved is not compu- 
tationally feasible. Thus approximate numerical and analytical solutions are useful. Such a 
solution for a static rotating torus was found by Krolik (2007). Making a number of assump- 
tions, he was able to obtain a semi-analytic model which showed that a rotating, static, and 
geometrically thick torus can be supported by infrared radiation pressure on dust grains. 

Dust opacity, is typically a ~ xlO times greater than the electron Thomson opacity, 
and thus the critical luminosity at which IR radiation becomes dynamically important is 
much smaller than the Eddington luminosity: L c ~ 10“ 2 — 10 -1 L e dd- If the temperature 
of the gas becomes larger than T vir>r then the radiation pressure prevails over gravity, and a 
model should include global plasma motions. 

In this paper we construct a model in which radiation pressure on dust grains not 
only supports the geometrical thickness of the torus but induces mass loss through infrared 
pressure driven winds. In our model, a ’’torus” is represented by an extended, dense and 
cold wind rather than by a static gravitationally bound torus. It is interesting that the 
physical conditions in such a. wind resemble those in red super-giant stars (except for the 
rotation) where radiation from a static ’’core” supports an extended, slowly outflowing en- 
velope (Bisnovatyi-Kogan & Dorodnitsyn 1999. 2001). In such stars, the outflowing wind 
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is driven by radiation pressure in the continuum, including significant contribution coming 
from radiation pressure on dust. 

In what follows, we numerically solve the equations of radiation hydrodynamics which 
describe the infrared-driven wind. Our solutions strongly support the concept of a dynamical 
torus: Compton-thick obscuration in which the structure is determined by infrared-driven 
flows of a dusty plasma. 

The plan of this paper is the following: we begin with basic assumptions underpinning 
our model in Section 2; the onset of an outflow is analyzed in Section 3; in Section 4 we 
derive the equations of radiation hydrodynamics describing the torus, and discuss appro- 
priate boundary conditions; our numerical method is outlined in Section 5; and results are 
presented in Section 6. The paper concludes with the discussion of major results, validity of 
approximations adopted and the relevance of our model to a physical picture of real AGN. 


2. Dusty torus supported by infrared pressure 

A spherically-symmetric distribution of fully ionized plasma around a. central mass can 
be gravitationally bound if the luminosity of the central object is L < L e dd, where L^d is 
the Eddington critical luminosity 


Ledd = 47rcGM — = 1.26 x 10 45 M 7j (1) 

where k t - 0.4 cm 2 g -1 is the Thomson opacity due to electron scattering, and M 7 = 
Mbh/ (10' Mq). 

The inner parts of an accretion disk around a black hole, where most of the accreting 
gas potential energy is dissipated, generate copious X-ray and UV radiation. Exposure 
of the outer region of an accretion disk to such radiation can have a profound effect on 
its structure and dynamics. In the following the dust opacity is denoted as k. In the 
UV the opacity of a single dust grain is significantly greater than Kg( v — 6 x 10 3 kt), 
adopting dust grain sizes 0.025 — 0.25 /on (Mathis et al. 1977), and dust grain density of 
2 — 3 gem -3 . Assuming perfect coupling between the dust and gas, and a dust to gas mass 
ratio, 50 - 100 the critical luminosity for the dust-plasma mixture becomes significantly less 
than the Eddington luminosity: 

^cdust - 5 x 10- 4 - 0.01 L e dd . (2) 

For instance, if a cold slab of plasma is exposed t o unattenuated X-ray and UV radiation, 



a significant part of such radiation will be absorbed and reprocessed into infrared in a thin 
” photospheric” layer of thickness, Sl/Ri pc ~ 1.3 x 10 ~ 3 n^ 1 (Hereafter, where appropriate 
we denote y x to represent a quantity y scaled in terms of 10 a: units of the same quantity, y). 

In the infrared, the Rosseland mean opacity, k of dust in the temperature range 10 2 — 10 3 
K is approximately 10—30 times larger than that of the electron Thomson opacities (Semenov 
et al. 2003). The dust opacity determines the critical luminosity, 

T 47t cG Al . , , 

L c — ~ (0.03 -0.1) L edd , (3) 

K 

If F = L/Ledd > 1, a spherically-symmetrical distribution of dust would be promptly 
blown away from approximately the dust condensation radius, ~ 0.3 — 1.5pc, for typical 
luminosities of 10 4a — 10 46 erg s _1 (Barvainis 1987; Phinney 1989). 

However, the presence of an equatorial accreting flow changes the picture. As a result 
of the reprocessing of external X-ray and UV radiation the incoming accretion flow (which 
otherwise would be geometrically thin) is pumped up with IR radiation and becomes geo- 
metrically very thick. For example, it has been shown that a thin disk (thin torus) eventually 
puffs up due to reprocessing of the hard X-rays in 10-100-keV range (Chang et al. 2007). 

As it becomes sufficiently fat, the torus intercepts significant fluxes of soft X-rays and 
UV. Reprocessing of this radiation to the infrared domain further pumps the torus interior 
with infrared photons, which become a major driving force in supporting the torus against 
collapsing back into a thin disk state. 

In a dusty plasma, the total pressure consists of that of an ideal gas, P g and that of 
radiation, n 

p = p s + n, (4) 

where 

P g = —pTZT, n = aT 4 /3, (5) 

f^m 

and 7Z — 8.31 • l(T(ergK“ 1 g^ 1 ) is the universal gas constant, a — 7.56 • 10 -15 (ergl<~~ 4 cm -3 ) 
is the radiation density constant, and y m is the mean molecular weight. Given the great 
variety of physical conditions in dusty molecular gas we set = 1 throughout this paper. 
For simplicity, we do not consider models which involve clumping, such as those of Krolik & 
Begelman (1988); Beckert & Duschl (2004); Honig & Beckert (2007), and assume continuous 
distributions of dust and gas. 

The relative importance of radiation pressure is described by the parameter 8 = P. P ~ 
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T 3 \ _1 

10 3 — + 1 j . At densities n — 10' cm' 3 . P g ~ II at T ~ 80K, and the ratio PJ II rapidly 

Tl'j J 

decreases at higher T, becoming 0.33 at 100K, and 0.04 at 200K. In the regime we are 
interested in, a few x 100 < T < 1000K, so that P g <C II. Thus, in this paper we neglect 
gradients of the gas pressure in the calculation of the equilibrium and dynamics of matter. 

The radiation energy density in the region we are concerned with is mostly determined 
by infrared radiation. On the other hand one can completely ignore the contribution from 
the mass-density of radiation, p r = aT 4 /c 2 ~ 8.4 x 10 -24 (gem -3 ) as it is much smaller than 
the mass-density of the gas, p ~ 8.35 x 10 -18 n 7 (gem -3 ). 

In the simplified model considered in this paper (Section 4) we relax the condition of 
vertical balance of the torus and treat it as a wind driven by the radiation pressure on 
dust. Self-gravity and dumpiness are ignored altogether. The equatorial inflow is implied 
but not calculated. It is also implied that this equatorial accretion inflow replenishes the 
gas lost in the outflow, but we do not attempt to model such a connection. We assume the 
flow is axially symmetric. One of the integrals conserved along the flux surface (i.e. such 
a surface which embraces a constant mass flux) is the specific angular momentum, l, (e.g. 
Beskin (2009)), and we assume the foot-points of the streamlines are located at the equator. 
Solving the momentum equation along z, we take into account only the v z dv z /dz component 
of the v ■ V v term of the equation of motion. The z-component of the radiation force is 
calculated from a 2D distribution of the radiation energy density, E{z, R ) which is obtained 
from the diffusion equation. The latter is solved numerically in 2D adopting the flux-limited 
diffusion approximation. 


3. The onset of the radiation driven wind 

Before embarking on numerical calculations (Section 4), it is instructive to consider a 
static model of a rotating torus. We approximate it by a spherically-symmetric distribution 
of plasma and radiation occupying a wedge of an opening angle # 0 , and extending along 
spherical radial coordinate r. Local thermodynamical equilibrium is assumed throughout 
the torus. Let us assume, (in this section only) that the radiation flux F which is given at 
the inner edge of such torus diffuses along r, as L/(Anr 2 ) ~ -Dr dT/dr , where L is the total 
luminosity, and Dy — 4acT 3 /(3«:p) is the diffusion coefficient. Limitations of such a model 
are obvious but for crude estimates we assume that there is no departure from spherical 
symmetry and L is conserved. This model resembles that of Bisnovatyi-Kogan & ZeLDovich 
(1968) who analyzed the onset of the stellar wind driven by high atmospheric opacity in the 
case of a non-rotating star. Here we extend their analysis by adding rotation. 
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The onset of the wind can be approximately derived by considering the radial balance 
equation and the equation for dT/dr, at fixed 9: 


1 dP = GM ( fi \ 

p dr r 2 \r 3 sm 2 9j' 

dT 4 3np L 

dr ac Airr 2 ’ 

where 9 is an angle measured from the vertical axis z\ R is the cylindrical radius, l — f 1R 2 is 
the specific angular momentum. fl(R) is the angular velocity which we assume to be constant 
on cylinders of constant R. Dividing (6) over (7), and formally integrating at a fixed 9 over 
r, and applying boundary conditions P = 0, T = 0 at the torus boundary, we obtain: 


P = 






( 8 ) 


where ly. = VrGM is the Keplerian specific angular momentum, and T c — L/L c , and 
E — aT 4 is the radiation energy density. To perform the integration in (8), we took into 
account that P is a single argument function of r, and assumed that / 2 /(/£) = a 2 = const , 
i.e. a radial model is specified by 9 and a. Taking into account that p = (P — aT 4 / 3) 
from (8), we obtain 


p, T3 = ^-{l-T c -a 2 ,sm 2 e). (9) 

Substituting (9) into (7), and integrating, we obtain 

T = To + ^ (1 - r c - a 2 / sin 2 0) (10) 

47 Z \r r 0 ) 

Where To — T(ro) at some fiducial r 0 . From equation (9) we conclude that the specific 
entropy of the radiation-dominated gas. S — 4/3 a T 4 / p is constant at constant 9. In the 
absence of rotation, equations (8), (9). and (10) are reduced to the corresponding equations 
of Bisnovatyi- Kogan & Zel’Dovich (1968). In order for (9) to be meaningful, the condition 
for a static torus follows: 


r c < 1 - o 2 / sin 2 9. (11) 

Alternatively, at fixed 9, an outflow begins if the condition (11) breaks down. From (10) 
follows a critical angle 9 C = arcsin(a:/v / l — F c ), such that, at 9 < 9 C a static configuration 



with a — const is not possible. At the equator, to be static such a torus should be sub 
Keplerian. 

Even if F c < 1 — a 2 / sin 2 6 from (10) and from the condition that T should not be finite 
at infinity we get another condition for the absence of an outflow 


r 0 <%(i-r c 


a 2 / sin 2 9). 


(12) 


where 


(13) 


is the gas virial temperature. If (12) is violated a wind will occur due to a combination of 
thermal and radiation driving. 


Using (9) in (10) and adopting a similar line of arguments one can deduce another useful 
condition for an outflow: 


T 0 >T vi , r (r 0 )ry\ 

where 


(14) 


T 

± V 




312 (— tMl ) 1/4 _ 987 ( 


njMf 


) 1/4 K, 


pc 


' pc 


(15) 


is the second definition of the virial temperature which replaces (13) in the case of II P g . 
In the radiatively-dominated torus, if condition (14) is fulfilled then an outflow driven by 
the pressure of the radiation flux, F — — Dj dT/dT begins. 


The relation for T v - IItI , (15) does not contain opacity. Notice that in the diffusion ap- 
proximation, the radiation force, g Ta d = kF/c ~ nXdE/dl ~ (1 /p)dE/dl, where A is the 
photon mean free path. A = 1 /(up). In the optically thick case dE/dl ~ E/L aT e 4 ff /(cT), 
where L ~ r is the size of the system. From balancing g rad and gravity GM/r 2 , the scaling 
(15) is obtained. 

In the free-streaming limit F = cE = crT^ B , where a — ac/A is the Stefan-Boltzmann 
constant, and here T e s is the temperature of the photosphere or of a layer where external 
radiation is converted to IR. Adopting the same line of arguments as in the optically thick 
case the following relation can be obtained: 
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r vir ,flx = - 292 (-i ^-) 1 / 4 K. (16) 

V ar K ) r pc K 10 

Note that the same result can be obtained assuming dE/dl ~ E/X, which is also valid in the 
vicinity of the conversion layer (i.e. in thermalization layer). In the optically thin case the 
radiation pressure is determined by the anisotropic radiation flux and in the optically thick 
case by the gradient of E, which is determined by the size of the system. The two effective 
temperatures are connected by the relation: 

d~virMx Tvir.n 

Td 

where t<± = Kpr is the optical depth parameter. 

The effective temperature of the conversion layer is found from oT F e dd 
is related to the total BH luminosity and Thomson opacity): 


(17) 

= cT 4 ff (here F 


T c 


off 


4 7 T 


GM 

K^cir 2 


1/4 


463 (£2^Ml)V4 Kj 


(18) 


pc 


where 7 0.5 is the fraction of the incident flux reemitted in the IR inside the torus. 


4. 2D + ID model 

Now we describe the ingredients of our numerical model for the radiation and gas flow. 
Consider equations describing stationary, slowly (v <C c) outflowing wind. The equation of 
motion and the continuity equation read: 

v • V v = G m - V$ 

V • {pv) = 0, 

where 

Gjr = —K, (21) 

c 

is the radiation force, and $ — —GM/(z 2 + R 2 ) 1 / 2 is the gravitational potential. For sim- 
plicity. w r e do not differentiate between the Rosseland and flux mean opacities (Mihalas & 
Mihalas 1984), and set dust opacity, k constant. We adopt a diffusion approximation which 
connects the infrared radiation flux, F with the infrared radiation energy density E: 


(19) 

(20) 
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F = — — X7E = —D VE. (22) 

3 Kp 

Notice, that (22) is in the form of the Fick’s diffusion law. The diffusion coefficient is 

D = cX v (23) 

In this paper we do not consider external heating by hard X-rays, and thus in the bulk 
of the flow we have: 

V • F = 0. (24) 

At small optical depths (i.e when p — » 0), the standard diffusion approximation breaks 

down: the mean free path, A — > oo, and D — » oo, and F —)• oo instead of F — >■ cE as 

it should be in a free-streaming limit. To take into account regions of r < 1 we adopt 
the flux-limited diffusion approximation (Alme & Wilson 1974; Minerbo 1978; Levermore & 
Pomraning 1981). In the flux-limited diffusion approximation A is replaced by A* = A A, 
where A is the flux limiter. The flux limiter we adopt is that of Levermore & Pomraning 
(1981): 


+ (25) 
where Alp = A \VE\/E. If r — » 0, then i? L p — > oo, and \F\ ~ cE. In the optically thick 
limit Alp — » 0 and A — > 1/3. 

Adopting cylindrical z, R coordinates, and assuming axial symmetry ( d/dcj) = 0) in the 
4> direction we numerically solve equation (24) in two dimensions. 

In z, R coordinates equation (19) takes the form: 

e z {v d z v + z P.2) + e R (R Q 2 k - R Q 2 ) = ~D VE, (26) 

where e : e R are coordinate unit vectors, and d x y = ||. We allow for only the v(z,R) z 
component of the velocity (hereafter v = v z ), i.e. an outflow is occurring along cylinders 
of constant R (see the end of Section 2 for discussion). The specific angular momentum, l 
must be conserved along the flux surfaces: l(R) = O (R) R 2 , and thus the angular velocity, P 
is constant on cylinders of constant R. The geometry of the flow is shown in Figure 1. The 
Keplerian angular velocity is found from 
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2 __ GM BU 
k ~ (R 2 + z 2 ) 3 / 2 ' 

Along the cylindrical flux surface, the continuity equation (20) reduces to simply 


(27) 


pv — (ir — const. (28) 

The amount of matter transported in ^-direction, p R is a function of the streamline, i.e. 
of R. It is convenient to rewrite the above equations in dimensionless units: x = R/Ro, 
z = z/Ro, E = E/Eq, p = p/po, v — v/vq (to simplify notation in the following we omit the 
tilde), where i?o, Eq, Po, and Vq are fiducial quantities. 



Fig. 1. — Illustration of the flow geometry and a coordinate system implied by the calcula,- 
tions (see the text for details). Not to scale. 

From (26), the momentum equation is cast in the following form: 

pv — — -A A x A 2 d z E - A 2 z nip, (29) 

az 

where we introduced dimensionless fij) = l/(x 2 + z 2 ) 3 ^ 2 , and the non-dimensional parame- 
ters, Ax, A 2 : 


A Eo 4 ggogo 

1 vi 0 n 2 oPo' vl ■ 

From (30) it follows that Ax = T|/7A r r . and A 2 = i^ 0 /ug, where v^o — 
r vir , r is found from (15). 


(30) 

(GM BIl /R 0 W 2 , and 
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In the current work we are not using the i?-component of (26) and cannot deduce the 
distribution of j in the moving wand self- consistently. That is because we are forcing matter 
to flow along cylinders of constant R and the radiation flux is pushing from just one side. 
The numerical solution which we obtain in Section 6 demonstrates that gravitation cannot 
balance centrifugal and radiation pressure forces and that departures from purely vertical 
motion should occur. 

The non-dimensional continuity equation reads: 


A* = P^i (31) 

where /i = Pr/(p 0 v 0 ), he. p(R 0 ) = 1. 

In order to solve equations (22), and (24) we need to find the distribution of density, 
p(z,x). This can be done solving equation of motion (29) and then using (31). It is conve- 
nient to convert equation (29) directly to the equation for p, making use of (31), and then 
numerically integrate this equation along the streamline. Thus, the equation for p reads: 

g = AAE + pz!®. (32) 

where d z E is known from the solution of the diffusion problem. Thus, equations (22), (24), 
(32) describe our problem. We emphasize that, although our treatment of the gas dynamics 
is quasi-one-dimensional (i.e. a flow along cylinders), our treatment of the radiation is fully 
two-dimensional. 


4.1. Boundary conditions and parameters governing the flow 

We solve equations (22), (24), (32) numerically in cylindrical coordinates z,R. In these 
coordinates, the computational domain has a rectangular shape with one side spanning from 
Rq to R\, and the other from z = 0 to z\. 

At the left boundary, we specify the distribution of energy density 


E(z, xq) = E x0 z e , (33) 

where E Xo = E(z — 0, Rq)/Eq. From the wind physics perspective, the case of smaller e 
mimics the situation when energy is deposited into the flow from the boundary over a longer 
region. 
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At the equatorial plane, at z = z 0 the flux is calculated from: 

F z (z 0 , R ) - -D dE/dz = aT* s , (34) 

where T e g is calculated from a ” photospheric” boundary condition, i.e. when T e g = T(z 0 , R) 
is obtained self-consistently when solving the 2D diffusion problem for E. 

At the upper boundary we apply a free-streaming boundary condition: \F\ ~ cE. We 
tried several implementations of the boundary conditions at the right boundary to find that 
the solution is not sensitive to their particular choice. However, it is reasonable to assume 
that the torus is close to being isothermal at larger R and not too large z, and thus we pick 
’’zero flux” boundary conditions at Ri. 

In order to obtain the distribution of p(z. Ri) on a particular flux surface, one needs 
to specify p(Ri). If we would have to match a stationary outflowing solution with a static 
solution in the accretion disk (i.e. vertical distribution of p) the situation would be equivalent 
to that described in Bisnovatyi-Kogan & Dorodnitsyn (1999): having at hand the vertical 
distribution of p in the accretion disk one would smoothly match it with the corresponding 
wind solution. This should be done at an arbitrary point z o, provided v z (zo) v s (zo), where 
v s is the sound speed, and from that matching the unique value of p w T ould follow. In our case 
we specify p(z — 0, R), and we must also specify p[R) (or v(z = 0, R)). At the equatorial 
plane we specify power law distributions for p and p: 

p(0,x)=x~ d , (35) 

We also choose that p scales as density at the equator, p(x) = x~ d , to provide v z (z = 0, x) = 
Vo is the same at all x in the equatorial plane. 


5. Solution: outline of the method 

Combining equations (22) and (24) we obtain the diffusion equation: 

V • F = — (D(VE) Z ) + ~(xD(VE) x ) = V r (E) + V J (E) = 0, (36) 

Equation (36) is solved numerically adopting an alternative direction implicit scheme (ADI), 
i.e. (Fletcher 1988; Fedorenko 1994). Here we outline the method while the details are left 
to Appendix B. 
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The computational domain {zi,Xi}, where i — 1, JV*, and j = 1 , Nj , spans from 0 to z\, 
and from Xo to x\ respectively. In our calculations, we adopt a 100 x 100 numerical grid 
which spans the 0 — 0.1 — 2 range in the z direction, and the x — 1 — 3 range in the x 
direction. We make use of a staggered grid: quantities E, p, and v are cell-centered, while 
D is face centered, (c.f. (Turner k Stone 2001)). In order to avoid approximation errors 
near the coordinate singularities when finite differencing in the curvilinear coordinates, we 
introduce volume elements dl% — xdx and dl\ = dz ( dl\ — dz is introduced for consistency), 
e.g., (Stone k Norman 1992). 

In order to solve equation (36), we introduce a pseudo time variable, t, and convert this 
equation into a time-dependent one: 


d t E-V 1 {E) -V 2 {E) = 0, (37) 

where d t y, V 1 , and V 2 schematically represent finite difference operators over t and along 
the alternative directions. 

In the ADI scheme, a single time step from t to t + St is made in the following manner: 
1) outer loop along 1 st coordinate (for example), with half time-step St* — St/2, with fully 
implicit scheme for the V l (E). Schematically, we have: <9 t (K* , E t j) — V 1 (E*) — V 2 {E) = 

0, where E* = (E*_ hj , Ef $t E* +Vj ), and E - (E hj ^ 1: E hJ , E lj+l | and E* = E(t + St*), 

1. e. applying a three-point stencil in a fully implicit numerical scheme for the update in 
1 direction. 2) Finally, iterating the outer loop in 2 nd direction: d t (Eij,E*,) — D 1 (E*) — 

V 2 (E) = 0, where E = Eij, E t j+i), and E = (£(_ij, E^j, E i+ ij), and obtaining 

E = E(t + 8 T *). 

Diffusion coefficients are taken at the ’’old” time, which has a tremendous benefit com- 
pared to dealing with linearized equations as would otherwise be necessary (in a fully implicit 
method). The fully implicit approach to the solution of a flux-limited diffusion problem was 
taken, for example by Hayes et al. (2006). 

As a consequence of a three-point finite differencing stencil implied by the diffusion 
operator in (36), the corresponding matrix equation for the updated E* and E involves a 
tri-diagonal matrix. We adopt a sweep method in order to solve the resultant tri-diagonal 
matrix equation via a tridiagonal matrix algorithm (Fedorenko 1994). 

A finite difference representation of the boundary conditions (BC) is derived in a way 
that preserves 2 d order accuracy of the numerical scheme. I 11 the ADI method, one can 
apply a combination of flux and temperature BC (Fletcher 1988). However the flux at the 
boundary should be parallel to one of the coordinate lines (Fedorenko 1994). In our model 
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the inner boundary is parallel to 2 and the flux should be normal to that boundary. Given 
our ignorance of the structure of the conversion layer at the inner boundary, we believe that 
is is slightly more physical to specify temperature BC instead of the flux one. Thus, for 
simplicity we specify the distribution of the effective temperature at the innermost cylinder, 
which marks the inner boundary of the computational domain. 

After the distribution of the radiation energy density, E(z, R) is obtained, the next 
approximation for p(z,R ) is found from (32). We solve this equation along cylinders, in a 
£ direction, adopting a 4 th order Runge-Kutta method (Press et al. 1992). The updated 
distribution of p is used to compute diffusion coefficients from (23) and again to solve (36), 
etc. The cycle is repeated until a stationary wind solution is found. 


6. Results of the numerical model 

It has been shown that the effective temperature of the conversion layer scales approx- 
imately as 463 (r 0 .5M7/7’p C ) 1 / 4 K, which is greater than ~ 312K and ~ 292K for 
a il/ 7 = 1 BH and no = 10 5 . It is reasonable to expect that no equilibrium is possible be- 
tween radiation pressure and vertical component of gravity and that a dynamic, outflowing 
atmosphere is a better description of what is going on. 

Our boundary conditions do not provide optimum acceleration as the incident radiation 
is normal to the flow at the boundary. It is the readjustment of the radiation flux inside the 
torus that produces a vertical gradient of E. Since in our simplified method we can calculate 
only the v z component of the velocity it is quite possible that taking into account the full 2D 
picture can increase terminal velocity (VE has its largest component approximately parallel 
to spherical r). The parameter n 0 scales the wind loading density. This density can be 
significantly smaller than the density at the equatorial part of the accretion disc. 

It is instructive to compare Vq ~ 210 ( Mi/R pc ) l > 2 kms -1 , with the sound velocity in the 
radiatively dominated plasma. Notice that in a radiation-dominated plasma v s ~ y'E/Sp, 
i.e. its value explicitly depends on both density and temperature. For relevant parameters 
we obtain: v s ~ w s ,rad — /ny) 1 ^ 2 kins" 1 . It is important that the wind launching speed 

is subsonic, and we choose vq = 0.1v s for all models. Note that since v s depends on p(z = z 0 ) , 
vq depends on it as well. 

We parametrize our models by the Thomson optical depth of the torus, tt = ktP dR 
calculated at the inclination 90° from the z-axis adopting the equatorial distribution of 
density (35) with d = 0.5 throughout all of the models. At the left, boundary we choose 
E x 0 — 1, and e — 0.1 in (33) and use T as a parameter instead of T e g which is calculated 
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from (18). 

The mass of the black hole is Mbh — 1 x 10' M®, and R 0 = lpc, and k = 10 are fixed 
for all models. 

Our results and various parameters of the models are summarized in Table 1. In the 
following we describe several characteristic models from the above set. 

We are not able to calculate models for F < 0.1 due to the intrinsic incapability of our 
method to treat low-velocity, decelerated flows. Such models require a full time-dependent 
multi-dimensional, radiation- hydrodynamics treatment. 

Models with the characteristic BH luminosity as low as 0.1L e dd produce a noticeable 
wind provided the optical depth is not too high. At larger optical depths the characteristic 
temperature a.t the equator is too low. As a result, we do not obtain an outflow solution for 
F = 0.1 and tt > 0.6. Increasing the optical depth from tt cs 0.2 to rx — 0.5 doubles the 
mass-loss rate to approximately 2M®yr~ 1 but also reduces the maximum velocity, u max by 
a factor of two. Most of the gas does not reach U esc forming a failed wind. With increasing 
tt the kinetic luminosity drops by an order of magnitude to L^ n ~ 1.2 • 10 3S (ergs~ 1 ), -which 
is an order of magnitude smaller than that obtained from simple estimates of the kinetic 
luminosity: Ly m ~ Mw T 2 nax /2 ~ 9.7 • 10 39 (ergs" 1 ). This is because only a fraction of the 
domain is occupied by the fast wind. 


Model 

r 

Ro 

r T 


n 0 

^max 

L 

rin 

Tbol 

M 

1 

0.1 

1 

0.17 

1 

• 10 5 

215 

4.22 

• 10 39 

1.24 • 

10 44 

1.23 

2 

0.1 

1 

0.34 

2 

■ 10 5 

153 

2.24 

■ 10 38 

1.24 • 

10 44 

1.74 

3 

0.1 

1 

0.51 

3 

• 10 5 

123 

3.59 

■ 10 38 

1.24 • 

10 44 

2.14 

4 

0.3 

1 

0.17 

1 

■ 10 5 

311 

1.63 

• 10 4 ° 

3.74 • 

10 44 

1.59 

5 

0.3 

1 

0.34 

2 

• 10 5 

217 

3.98 
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3.74 ■ 

10 44 

2.25 

6 

0.3 

1 

0.51 

3 

• 10 5 

251 

1.33 

• 10 39 

3.74 ■ 

o 

T 1 

2.76 

7 

0.3 

1 

0.85 

5 

■ 10 5 

129 

5.34 

• 10 38 

3.74 • 

10 44 

3.56 

8 

0.5 

1 

0.17 

1 

• 10 5 

445 

4.65 

1 — 1 
o 

o 

6.24 ■ 

o 

r— 1 

2.05 

9 

0.5 

1 

0.51 

3 

• 10 5 

357 

5.9 • 

t— 1 

o 

O 

6.24 • 

10 44 

4.1 

10 

0.5 

1 

1.48 

5 

• 10 5 
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1.2 ■ 

10 4 ° 

6.24 ■ 

10 44 

5.29 

11 

0.5 

1 

2.37 

8 

• 10 5 
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2.76 

■ 10 39 

6.24- 

10 44 

6.69 

12 

0.8 

1 

1.48 

5 

• 10 5 
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4.2 • 

10 4 ° 

9.99 • 

10 44 

6.7 

13 

0.8 

1 

2.37 

8 

• 10 5 
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1.04 

10 4 ° 

9.99- 

10 44 

8.47 

14 

0.8 

1 

2.97 

1 

■ 10 6 

162 

5.56 

10 39 

9.99- 

10 44 

9.47 

15 

0.8 

1 

5.94 

2 

• 10 6 

104 

1.19 

10 39 

9.99 • 

10 44 

13.4 


Table 1. Models characterized by different initial parameters: F, i?o(pc), tt, charac- 
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teristic density n 0 (cm 3 ), and resulting kinetic and bolometric luminosities, Tkin(ergs 1 ), 
Lboi(ergs _1 ), and mass-loss rates M(M 0 yr -1 ). 

The density and radiation energy density for Model 1 are shown in Figure 2, and Figure 
3 shows the surface plot of the velocity, v/U esc , where U esc is the local escape velocity. Recall 
the distributions of E and p at the appropriate boundaries (33). The most appropriate 
conditions for the acceleration of the wind happen in the middle of the domain, where the 
radiation field has strong gradients, but density is lower than that at the left boundary. The 
maximum value of the effective temperature T 0 = T(R 0 ) = 407K which rapidly declines at 
larger spherical radii, r. Maximum velocity attained by the wind is 4.7 Ma, where Ma is the 
Mach number. Further increasing tt results in a drop of the velocity: Most of the wind has 
velocity smaller than escape velocity, however the ’’mass-loss rate” of such a failed wind is 
noticeably larger (c.f. Table 1). 

From Figure 2 (right panel) one can see a significant drop of radiation energy density, 
E within the distance, Sx ~ 1.5 from the left boundary and from Figure 3 we identify the 
region x ~ 1.5 — 2.5 and where the most of the fast wind is blowing. 

As was discussed in Section 4.1. in the approach taken in this paper, the mass-loss rate 
from the two sides of the disc, M is a mere consequence of the adopted boundary conditions. 
From the continuity equation (31) p ~ p/v z , and thus an increase of the velocity in region I is 
compensated by the reduced density in accord with what is observed in Figure 2 (left). The 
enhanced density region in Figure 2 (left panel) corresponds to a low velocity, high density 
and quasi-isothermal region of the torus. In the following we denote the high velocity part 
as region I, the higher density, narrow transition region as region II and the high-density 
region located at larger radii as region III. 

Models 8-11 have F = 0.5. These are optically thin, marginally optically thick and 
Compton thick models. Increasing tt from ~ 0.2 to ~ 1.5 results in increasing M from 
~ 2M 0 yr“ 1 to ~ 5.3M e yr 1 . The maximum velocity drops from u max ~ 445 km s 1 for 
Model 8. to ~ 203 kms -1 for Model 10. The color intensity plot of p(z,x ) and E(z,x) are 
shown in Figure 4. One can see the fast wind occupies approximately 40% in the radial 
extent, and that in the wind region the density is markedly lower then in the outer parts. At 
larger R there is no significant outflow. Notice the locus of a sharp rise of the density which 
marks a barrier between the region of fast flow r and the almost quasi-static torus. The wind 
is supersonic, for example, the maximum Mach number for Model 9 is 6. 

Further increase of tj to 2.37 engages more matter into the low velocity wind. The color 
intensity plots of p(z.x) and E(z,x) are shown in Figure 5. The wind in Model 11 does not 
reach the local escape velocity: the maximum velocity is 0 .6 U esc (z . R) . Large amounts of 
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Fig. 2. — Model 6: Color-intensity plots of the dimensionless density, p (left) and dimen- 
sionless infrared radiation energy density, E (right). Axes: distance in parsecs. 



Fig. 3. — Model 6: Velocity surface plot. Axes: vertical: velocity v z /U esc , where U esc is 
the local escape velocity: horizontal: R : distance from the BH in parsecs; 2 : distance from 
equatorial plane in parsecs; 
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Fig, 4- - Model 10: Left: velocity surface plot, where v z /U e sc , and U esc is the local escape 
velocity. Right: density surface plot. Horizontal: R: distance from the BH in parsecs; z: 
distance from equatorial plane in parsecs; 

gas, 6.69M 0 yr _1 participate in a low velocity ~ lOOkms -1 flow, leaving the computational 
domain in the form of a failed wind. The high-density region III is clearly seen in Figure 5 
(left panel). 

Models 12 -15 in Table 1 represent a BH shining close to L e M- They have F = 0.8; 
These models are Compton-thick: rj — 1.5 — 6. Results for Model 13 are shown in Figure 6. 
and for Model 15 in Figure 7. From Figure 5 and Figure 7 one can see that within a high 
density region there is a higher density '’core", which is most pronounced in Model 15. 

The high velocity wind in Model 15 occupies a narrow wedge-like region close to the 
left boundary, but even there, the maximum velocity, v max ~ 0.5U esc (z, R). Low velocity 
of the wind translates into a ratio L\\ n / L\- JO \ dropping to ~ 8 ■ 10~ 7 . Most of the flow is 
mildly supersonic, Ma~ 1.5 — 2 with velocity below the escape velocity. Here AI is better 
interpreted not as a mass-loss rate but as a parameter describing how much gas is involved 
in large-scale motions; its value approaches 13.5M Q yr _1 . 

Most of our models demonstrate a clear separation of the torus into lower and higher 
density parts. From arguments of Section 2 one can expect that some sort of a transition 
region should exist between a quasi-isothermal torus ’’core" and an infrared-driven part 
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Fig. 5. — Model 11. Color-intensity plots of the dimensionless density, p (left) and dimen- 
sionless infrared radiation energy density. E (right). Axes: distance in parsecs. 

located closer to the source of UV and X-ray radiation. The existence of the over-dense region 
in our 2D gas distributions further supports this idea. One can argue that an interesting high 
density region observed in Figure 6.7 can be a sign of a quasi-static/stationary core. Such a 
region would be a likely place for large scale meridional motions which wrap a quasi-static 
region. The final answer can only be provided by a fully 2.5D time-dependent radiation- 
hydrodynamical simulations. 


7. Discussion and Conclusions 

We have studied a model of an AGN torus in which obscuration is provided by a 
radiation-driven wind rather than by a static distribution of gas. This can occur if the 
UV and X-ray radiation generated in the inner parts of an accretion disk is reprocessed into 
the infrared (IR) in the cold, dusty environment at approximately lpc from a supermassive 
black hole. We have shown that due to high dust opacity the pressure of such IR radiation 
has a profound effect on the torus dynamics and structure. 

Semi-analytic models of a static rotating torus which is supported by infrared pressure 
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Fig, 6. Model 13: Left: velocity surface plot. Right: density surface plot. Axes: vertical: 
left: velocity v z /U esc , where U esc is the local escape velocity; right: density: horizontal: z: 
distance from equatorial plane in parsecs; R: distance from the BH in parsecs; 



R 


Fig. 7. - - Model 15. Color-intensity plots of the dimensionless density. Axes: distance in 
parsecs. 
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on dust grains have been developed by Krolik (2007). The results confirmed that the torus's 
thickness can be entirely supported by IR radiation, and also raised new questions. One of 
the important ones is how to construct a model in which the plasma avoids being blown away 
(Notice, that critical luminosity L(;J ust <0.1 L e dd.e); without fine tuning of the parameters? 
In this work we relax the assumption of the static torus and suggest a model which takes 
into account plasma motion. 

We adopt several approximations and simplifying assumptions. Only constant dust 
opacity was taken into account, despite the fact that close to the wind (torus) surface a 
significant portion of the dust should sublime due to X-ray heating. In our model we assumed 
that the wind possesses only a vertical component of the velocity by forcing it to move along 
cylindrical surfaces. In reality we expect the wind to become more radial at large r. 

We neglect self-gravitation of the torus despite geometrical and column density argu- 
ments which favor torus masses of 10 4 — 10 5 M e . It is likely that the self-gravitating instability 
may operate inside such a torus forming an interacting system of molecular-dusty self- grav- 
itating clouds (Krolik & Begelman 1988: Beckert & Duschl 2004). If there is enough column 
density and the torus is already geometrically thick it will inevitably intercept and convert 
UV and soft X-ray radiation into IR, providing vertical support and possibly suppressing 
the self-gravitation instability (Thompson et al. 2005; Honig &: Beckert 2007). Thus the 
optical thickness tt of the torus plays important role as an optically thin self-gravitating 
torus would likely collapse into a thin disk with subsequent star formation (Toomre 1964). 

In our simplified model, the torus is described by equations of continuous radiation 
hydrodynamics. Even such an oversimplified approach required a complicated numerical 
treatment. The most important part is that in order to obtain the distribution of radia- 
tion energy density, E we numerically solve a 2D diffusion equation adopting a flux-limited 
diffusion approximation. 

In our method we solve a simplified system of equations of radiation hydrodynamics, 
assuming a stationary outflowing wind driven by gradients of IR radiation pressure. Our 
method is not free from serious limitations: we cannot follow ’’marginal” situations, such 
as a slowly outflowing wind with deceleration. For example, if somewhere in our 2D com- 
putational domain such a situation happens, the calculation must stop. This happens, for 
example, if the BH luminosity is too low, < 0.1L ec jd, or the density is too high, tt > 6 (see 
Table 1). 

In most of our simulations we find three characteristic regions: In region I conversion of 
external UV and soft X-rays into IR provide ample radiation pressure not only to support 
the torus vertically but to initiate a rigorous outflow; In region I, the radiation pressure is 
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strong enough to accelerate plasma to velocities 300 — 400 kins -1 (for tt cs 0.5, T > 0.5): At 
larger R a narrow region II is located where the density rises and the wind is either failed 
(i.e. first accelerated and then decelerated) or decelerated. Region II acts as a barrier 
separating the dynamical part from the quasi-static one. The latter we call region III and 
velocities and densities there are small. Region III is quasi-isothermal, although the vertical 
gradient of the radiation pressure is large enough to support its geometrical thickness. 

In a real torus, the global flow pattern should be very complex: soft X-rays heat the torus 
surface where cooling cannot compensate for the radiation heating, the temperature rises 
sharply and outer torus layers start evaporating. Numerical simulations (Dorodnitsyn et al. 
2008), show the formation of a wind with temperature, T w ~ 10 4 — 10 'K which evacuates 
10~ 3 — 0.1 M© yr _1 from the torus. In the bulk of the warm absorber flow radiation pressure 
plays almost no dynamical role. Multiple phases of the cold and hot gases may co-exist in the 
outflow (Krolik et al. 1981). In addition to the evaporated gas, UV-line-driven winds (Proga 
et al. 2000) wdiich are stripped from accretion disk at much smaller radii can also contribute 
to filling the funnel of the torus. The incident UV and soft X-ray flux is attenuated in this 
gas. In the bulk of the torus the gas pressure is much smaller than the pressure of the infrared 
radiation which is the major force which keeps the torus geometrically thick. The gradients 
of gas pressure become important in the narrow evaporative layer where temperature jumps 
from the cold inner values to the values corresponding to the temperatures of the warm 
absorber gas, T w . 

The very high opacity of the cold dusty plasma completely stops UV radiation some- 
where further into the torus, within the narrow' layer of the UV photosphere. In our current 
work we identified such a UV photosphere with the inner torus boundary. The radiation 
input was prescribed assuming the distribution of the effective temperature at this boundary. 

Further from the photosphere rotation plays an important role in shaping the density 
and IR optical depth contours. The locally super-critical IR flux creates a radiation-driven 
outflow. The conditions in such a flow resemble those in winds of supergiants (Laniers & 
Cassinelli 1999) or evolved massive stars (Bisnovatyi-Kogan & Dorodnitsyn 1999). Hard 
X-rays, with energies E > 10 keV, penetrate much deeper into the torus body creating 
significant local deposition of energy through Compton scattering (Chang et al. 2007; Shi 
& Krolik 2008). Some contribution to the radiation field can also come from star formation 
taking place within the torus or the obscuring flow' (Wada & Norman 2002). We will study 
the influence of these important effects in a future paper. 

Close to the torus boundary the infrared flux, Fir ~ — VT, propagates approximately 
along the inside normal to the surface. The curvature of the photosphere will significantly 
influence the distribution of temperature in a thin (of the order of a few mean free paths) 



- 24 - 


thermalization layer. Deeper into the torus, the rotation, the gravitation force, and plasma 
motion (through the continuity equation) determine the distribution of p. At higher heights, 
2 : density tends to be lower, and VT is more parallel to Vtir, where tir is an optical depth. 
As a result the infrared radiation diffuses in a direction ~ — Vtjr. Further into the torus 
the radiation tends to make it isothermal. Inside these regions there still exists a significant 
component of the radiation pressure in the 2 direction, but the quasi-static approximation is 
applicable and the torus is described by models such as those of Krolik (2007); Shi & Krolik 
(2008). In the intermediate region convective transport of energy may be of importance. 

The torus loses mass with an average rate of M ~ 1 — 10 M© yr~ 4 with negligible kinetic 
luminosities, L kin which are <C 1 % Lboi depending on various model parameters. This leaves 
two possibilities: If the gas escapes from the system then the torus will be depleted within 
10 4 — lO’yr which brings an important connection of the IR-driven obscuration with the 
AGN feedback problem. Taking into account that radiation-driven flow's are believed to be 
important for the AGN feedback, see e.g. (Begelman 2004; Fabian 2010) our results may 
further favor these ideas. 

Yet another possibility is that a considerable part of matter does not leave the torus’s 
potential well, instead forming global vortex-type motions. The impossibility of balancing 
in one static picture radiation, gravitation and a centrifugal forces is a well known cause of 
meridional flows in rotating radiative stars (Tassoul 1978; Kippenhahn & Weigert 1994). In 
a thin accretion disk such imbalance leads to the mass outflow from the disk at a luminosity 
considerably smaller than the critical one (Bisnovatyi-Kogan & Blinnikov 1977). 

A high density component in the meridional cut of the velocity distribution is often 
present in purely hydrodynamical simulations of accretion flows and winds. For example, 
such purely hydrodynamical 2.5D simulations of Dorodnitsyn et al. (2008), (Figure 5) reveal 
the presence of meridional- like, returned current. Notice that such a region is also present 
in Figure 7 of our simulations as well. In these works energy transport is performed by 
advection. As was shown here, it is entirely plausible that II P g in the bulk of the 
obscuring flow. Thus, inclusion of the infrared radiation pressure and radiative diffusion and 
advection of the radiative energy density into a time-dependent hydrodynamical framework 
should demonstrate wdiether the bulk of the torus is quasi-static and IR supported, or in a 
form of a dusty IR-driven flow. 
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Appendix: solution of the diffusion equation 

The diffusion coefficient is defined on the augmented grid which is shifted from the i,j 
grid by h x along the i axes to the left and by h y down along the j axes. Thus, with regard 
to i,j cell we have: D\ } is located at the left i boundary; Dj +lj at the right i boundary; 
D'lj at the downside j boundary; D 2 J+1 at the upper j boundary. 

When making a time step r, in the ADI scheme we first perform the inner sweep along 
i and the outer along j and then alternate i and j inner and outer sweeps. First we calculate 
E* = E(t + t*), and then E = E(t + r). If the inner sweep is along the I th index, then the 
finite difference equation to be solved at i,j reads 
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This preprint was prepared with the A AS MpK macros vo.2. 
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where r* = r/2. Then, the inner sweep is made the j th index, and the outer along j. The 
finite difference equation to be solved at i,j reads 
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A sweep method is adopted to solve the resultant tri-diagonal matrix equation via a 
tridiagonal matrix algorithm (Fedorenko 1994). 

A finite difference representation of the boundary conditions (BC) should preserve 2 d 
order accuracy of the numerical scheme. For example, if the zero flux BC are given at the 
inner i boundary, we write 


{Eis+ij - E is -.i t j) / (2h x ) = 0, (42) 

where is is the first index along i, and is — 1 is the index of the ghost zone. The idea is 
to express E is ^ij at the ghost zone from the relation for the BC such as (42), and then 
to substitute the result into equation (38) written for the is zone. The resultant equation 
couples only is and is + 1 indices: 


where 



(43) 


Slj = t | (Di J+ 1 (E iSij+1 - E is j) - Dl.{E iStj - Euj-x)) + E is j-f. (44) 

n y r 

Finite difference representations of the boundary conditions at other boundaries are derived 
in a similar fashion. 



